Multiply robust estimation of marginal structural models in observational studies subject to covariate-driven observations

ABSTRACT Electronic health records and other sources of observational data are increasingly used for drawing causal inferences. The estimation of a causal effect using these data not meant for research purposes is subject to confounding and irregularly-spaced covariate-driven observation times affecting the inference. A doubly-weighted estimator accounting for these features has previously been proposed that relies on the correct specification of two nuisance models used for the weights. In this work, we propose a novel consistent multiply robust estimator and demonstrate analytically and in comprehensive simulation studies that it is more flexible and more efficient than the only alternative estimator proposed for the same setting. It is further applied to data from the Add Health study in the United States to estimate the causal effect of therapy counseling on alcohol consumption in American adolescents.


INTRODUCTION
The study of causes and effects is an es s e n tial compone n t of le arning he althcare system s ( Krumho lz, 2014 ) .Estim ate d causal effects, in stead of as s oci ation s, should be us ed to inform tr ea tme n t de cisions .This m a n us cript propos e s a nove l , mult iply robus t efficie n t es tim ator for the m argin al causal effe ct of a tr ea tme n t that may vary in time on a longitudinal outcome using o bs erv ational d ata.Randomize d c ontrolle d trials are the gold s ta nda rd for causal infe re nce.Ra ndomization of the treatme n t options makes patie n ts more comparable in terms of their baseline cha racte ris tics.Ra ndomized s tudies also ofte n h av e clear protocols for the timing of patie n ts' visits ( ie, o bs erv at ion t imes ) a t which pa tie n t health s ta tus is measur ed.It is not alw ays pos sib le to conduct a r andomiz e d c ontrolle d trial designed to answer a spec i fic causal question and r esear chers often turn to o bs erv ational d ata ( Bl ack, 1996 ) .We herein focus on the pa rticula r fea tur e s of e lectr onic health r e c ords ( EHRs ) data.
While EHRs ar e incr easingly av ail ab le for analysis, they are not c olle cte d for r esear ch purposes.The tr ea tme n ts meas ure d in EHRs are not r andomiz e d to patients .This can lead to spurious as s oci ation s in the data called confo un ding ( Gree nla nd a nd Morge ns te rn, 2001 ) .Thes e d a ta ar e also measur ed irr egularly acr oss pa tie n ts.Each patie n t follows their pa t tern in how they ac c es s care, als o like ly to de pe nd on their cha racte ris tics.For example, as in B ůžková and Lumley ( 2005 ) , suppose an o bs erv ational study in which we are int erest ed in the causal effect of air pollution on forc e d expiratory v o lume ( FEV ) .We as sume the effect of air pollution on FEV is further me diate d b y as thma, and that both air pollution and asthm a affe ct the ch anc e for FEV to be meas ure d.Or, as in our application in this ma n uscript, s uppose w e aim to estim ate the m argin al causal effe ct of thera p y on the average number of alco ho lic beverages consumed in Ame rica n ado les ce n ts follow e d irre gularly ov e r time, whe re obs erv ation depends on ado les c ents' ch aracteristics .Statistically, this cr ea t es a long-t e rm depe nde nce s tructure betw e en the outcome and the visit proces s es that can result in bi as ed e stimators of causal or associational pa ra mete rs ( se e e g, Lin and Ying 2001 ;McCulloch et al. 2016 and most re c ently Coulombe et al. 2021Coulombe et al. , 2022 ; Ya ng 2022 a nd P ulle n aye gum et al. 2023 in a c ontext of causal infe re nce ) .Whe n aiming for a causal effect, this bias can be due to confounding by the visit process or, if the visit ind icator s act as coll ider s ( ie, are affe cte d by the tr ea tme n t prescribed and the study out come ) , t o collide r-s tra tifica tion bias ( Gree nla nd, 2003 ) .In the examples liste d abov e, the causal re lation ship be tw e en air pollution a nd FEV measure me n ts or the ra p y a nd alco ho l con sumption is also likely c onfounde d.
Under a s e t of caus al and modeling as sumption s, caus al effects ca n be infe rred b y es t imat ing the pa ra mete rs of a m argin al structural model ( MSM ) fitted on the da ta fr om a pseudo-population that is free of confounding a nd othe r types of spurious ass oci ation s, such as collide r-s tra tifica tion bias ( Robins et al., 2000 ) .Previous work has tackled this in s e ttings with cov ari atedriven o bs erv at ion t imes and confoundin g, lea din g to the flexib le invers e pro bability of tr ea tme n t a nd monitoring wei gh ted ( FI PTM ) est imator ( Coulombe et al., 2021 ) .For observation times occurring not at ra ndom, P ulle nayegum et al. ( 2023) rece n tly proposed a nothe r es tima tor tha t uses random effects to model the remaining depe nde nc e betw e en the observation and the outc ome proc es s es.Both me thods above may suffer similar iss ues: They w ere not dev elope d to be the mos t efficie n t es timators in their se mipa ra metric class, a nd they are not doubly r obust, but ra ther r e ly on mode l as sumption s for the tr ea tme n t and the o bs erv at ion t imes.Yang et al. ( 2020 ) and Rytgaard et al. ( 2023 ) als o propos ed s e mipa ra metric a pproaches for the timespec i fic in te rve n tion effects.Their a pproaches we r e pr oposed for the study of survival outcomes as opposed to con tin uous outc omes .Yang ( 2022 ) and Rytgaard et al. ( 2022 ) proposed ge ne ral se mipa ra metric fra me works for e st imat ing in te rve n tion-spec i fic mean outc omes .Most of these approaches allow several longitudin al proc es s es in the est imat ion ( exposure , outcome , and cov ari at es ) t o be measured sporadically by jointly modeling all their o bs erv a tion pr oces s es.They are highly flexible but the estima nds a nd es t imat ion approaches used by these authors are diffe re n t tha n our s.They foc us on the mean outcome d iffe re nce over a pre-spec i fied period of time, under a spec i fic tr ea tment r egime.We her ein pr opose a s trai gh tforwa rd, es tim ating e quation approach for the average tr ea tme n t effect es timated with r epea te d meas ure me n ts for which R code, av ail ab le with the ma n uscript, is s trai gh tforwa rd to imple me n t.We focus on the situation when the o bs erv at ion t imes occur "at random" ( as oppose d to c ompletely at random or not at random ) .The FIPTM es timator ca n be used in that s e t ting, but it r e lie s on the c orre ct spec i fication of the treatme n t a nd outcome o bs erv ation models as a function of patie n t cha racte ris tics, a nd ca n be seve rely bi as ed when one or both models are not c orre ctly spe cifie d.Se c ondly, the FIPTM could be made more efficie n t b y de riving the influenc e curv e for the causal effect of in te res t ( Tsiatis, 2006 ) .To addre ss the se issue s, we propose the first multip ly ro bust estimator for the causal m argin al effe ct of a bin ary tr ea tment on a longitudin al c on tin uous outc ome, th at ac c ounts for c onfounding and irre gular c ovar iate-dr ive n obse rvat ion t imes of the outcome sim ulta neously.The notation, es tima nd, caus al as sumption s, and propose d estim ator are prese n ted in Sect ion 2 .Simulat ion studies c ov ering sev e ral diffe re n t sce na rios of data ge ne rating mechanism ( DGM ) ar e pr esented in Section 3 .Our me thodo lo gy is app lied to d a ta fr om the Ad d He alth study in Section 4 and we conclude in Section 5 .

METHODS 2.1 N ot ation
We ass ume w orking with a ra ndom sa mple of size n from a la rge r population, i denotes the patie n t index a nd t ∈ [0 , τ ] is the time with τ a maximum cen s oring time in the co hort.Le t A i (t ) represe n t the bina ry tr ea tme n t t aking value s in { 0 , 1 } and Y i (t ) be the con tin uous outcome for patie n t i at time t.We de note v e ctors and m atric es in bold.The type of DGM we focus on is prese n ted in the left panel of Figure 1 , in which K i (t ) are potential confounders for the tr ea tment-outcome r ela t ionship ( Pearl , 2009 ) , M i (t ) are potential med iator s for the treatme n t effect on the outcome, and P i (t ) contains pure pred ictor s of the outcome th at c ould also affe ct the o bs erv a tion of a pa tie n t outcome Y i (t ) .
The s e t P i (t ) i s di stingui she d from K i (t ) as it c ould c ontain visit pred ictor s ge ne rated a nd meas ure d after the tr ea tme n t.Only the outc ome proc es s is as s ume d to be meas ure d sporadically ( e g, the wei gh t is meas ure d irre gularly ac c ording to patie n t cha racte ris tics such as a change in medication ) .All the other v ari ab les ne c es s ary to the est imat ion of the marginal causal effect of tr ea tme n t a re ass ume d to be av ail ab le at all times during follo w -up; in EHRs, the drugs and comor bid ities a re ofte n re c orde d anytime there is a new diagnosis or a new prescription, so it is often a reas onab le as sumption .Le t N i (t ) be a c ounting proc es s for o bs erv at ion t imes of the outc ome Y i (t ) betw e e n times 0 a nd t for individual i .The indicator d N i (t ) equals to 1 when there is an o bs erv ation of the outcome Y i (t ) , and 0 otherwise.The set V i (t ) includes all the v ari ab les causing o bs erv at ion t imes ( ie, causing dN i (t) ) and also includes all the confounders of the tr ea tme n toutcome r ela tionship.We m us t include K i (t ) in the s e t V i (t ) for the proposed estimat or t o be consis te n t.In Fi gure 1 , we have Patie n ts a re allow e d to h av e diffe re n t ce n s oring times C i .Denote ξ i (t ) = 1 { C i > t} an indicator of patie n t i s ti l l in the study at time t.We ass ume th at c en s oring times ar e non-informa tive, a n assumption de noted b y Y i (t ) ⊥ C i | A i (t ) ( thi s i s di s cus s ed more in Section 2.5 ) .

Assumption Definition
Outc ome c onsis te ncy re c eiv e d treatme n t option A i (t ) = 1 or A i (t ) = 0 , respe ctiv ely.
The tr ea tme n t m ay be fixe d at baseline or it c ould vary in time.In wh at follows, w e define A i (t ) to be the tr ea tme n t give n at time t.
T he causal mar ginal effect of the binary tr ea tment on a continuous outcome is defined as Our in te re st lie s in a cross -se ction al effe ct β 1 that does not vary in time, which can be estim ate d using an MSM with which w e ass ume a cons ta n t effect ( see the Dis cus sion s ection in which this assumption is d isc ussed ) .
Suppose a certain time dis cre t izat ion for which the re ca n be only one jump in the c ounting proc ess N(t ) ( for ins ta nce, daily visits at the doctor's office ) .If we had ac c ess to all pote n tial outcomes under both tr ea tments and at each time t for the time gra n ula rity chose n above, for a ra ndom sa mple of pa rticipa n ts of size n , we could estimate β 1 using s amp le mean s.On the other ha nd, b y conducting a r andomiz ed contro lled tri al and randomly alloca ting pa tie n ts to one of the two tr ea tme n t options, a nd obse rving patie n ts at prespec i fied visit times, then patie n ts allocat ed t o treatme n t 1 a nd tr ea tme n t 0 should not differ before re c eiving the treatme n T could be used to estimate β 1 without a ccountin g for confoundin g.In o bs erv ational da ta fr om EHRs, unfortuna t ely, we t end t o o bs e rve the pote ntial outcomes Y 1 i (t ) in those who had gr ea te r cha nces of being tr ea ted with A i (t ) = 1 , a nd the pote n tial outcomes Y 0 i (t ) in those who had gr ea te r cha nces of being tr ea ted with A i (t ) = 0 ( as a con s eque nce, In addition, the pote n tial outcomes for an individual i are only observ e d at times t when dN i (t) = 1 , which may depend on cov ari ates.We do not have access to all pote n tial outcomes a nd require causal assumptions to equate the es tima nd to functions of Y i (t ) .

Cau s a l ass u mpt ion s
Five causal assumptions ar e r equir ed for consis te n t es t imat ion of the causal m argin al effe ct of tr ea tme n t ( Table 1 corresponding to as sumption s 1-3b be low ) .Mode ling as sumption s for the MSM a nd the n uisa nce models are also r equir ed, the la t te r a re d isc ussed in Section 2.5 .These five assumptions m us t hold, both for the former FIPTM es timator a nd for the nov el propose d estim at or t o be consist ent, except for the inclusion of K i (t ) in the s e t V i (t ) that is only required for the proposed estimator.The conditional exchangeability can be re c ov ere d by breaking the spurious as s oci ation s due to the tr ea tme n t a nd o bs erv ation me ch anism s vi a invers e wei gh ts ( ma rginal a pproa ch ) , conditionin g on the s e ts A i (t ) and V i (t ) ( which includes K i (t ) ) in a r egr e ssion mode l for the outcome and using methods such as g-computation ( Rob ins, 1986 ) ( s ta nda r diza tion appr oach ) , or as we propose, using both approache s simult aneously to obtain a robust estimator.We review the previously proposed FIPTM estimator.

Pr evious estim ator
Using the m argin al approach c orresponds to using the FIPTM estima tor pr oposed in Coulombe e t al. ( 2021 ) .It con sists of a doubly-wei gh ted leas t squa re s e stim ator th at inc orporates invers e pro bability of tr ea tme n t ( IPT ) wei gh ts ( Horvitz a nd Thomps on, 1952 ; Ros enb a um and Rubin, 1983 ) and inverse intensity of visit ( IIV ) weights ( Lin et al., 2004 ) .The IPT wei gh ts are functions of the confounders K i (t ) and the IIV weights are functions of the visit pred ictor s V i (t ) .The estimator is consiste n t for β 1 when both weights are c orre ctly spe c i fied.A pa ra metric model can be used for where ; ψ } , the propen sity s core, is the pr obability of r e c eiving the tr ea tme n t 1 as a function of predictors K i (t ) and parameters ψ ( Rosenb a um and Rubin, 1983 ) .A log i stic r egr ession can be used to compute an estim ate d propensity score.The IIV wei gh ts ca n be obtained b y modeling the mean visit indicator as a function of cov ari ates V i (t ) .Since the vi sit indicator i s bina ry a nd visits a r e r ecurr e n t, one ca n use a model for r ecurr e n t visits such as the Ande rse n a nd Gi l l ( 1982 ) mode l ( which corre sponds to a proportional rate model ) or a log i stic r egr e ssion mode l.Both mode ls re ly on re lative ly similar as sumption s for the mean visit indicator when using the same s e t of cov ari a tes, but the pr oportional ra te mode l mode ls the rate and the log i stic r egr ession, the pr oba bility of visit.T hey lead to simila r es timates of the rate and probability of visit when the visits a re ra re s uch th at only one visit occ ur s over a time unit ( eg, a d ay ) , s ee e.g., Papoulis and Pi l lai ( 2002 ) .The pr oportional ra te means correctly spec i fied and X means no r equir e me n t.

Scenario
model for the visits as a function of V i (t ) is given by The baseline rate of o bs erv ation λ 0 (t ) in ( 2) consists of the visit rate when all v ari ab les V i (t ) are s e t to their refe re nc e lev el.With the FI PTM est ima tor, the baseline ra te can be r emov e d from the IIV wei gh ts without affecting the marginal effect of tr ea tme n t estim ate sinc e this w ould sti l l m ake the w ei gh ts in ( 3 ) proportional to the in te nsity of being observ e d as a function of V i (t ) .
The IIV wei gh ts ca n also be s tab ili zed, in which cas e the bas eline rate canc els autom atically in the wei gh ts a nd ne e d not be estim ate d ( B ůžková and Lumley, 2009 ) .This leads to the follo wing in te nsity of visit wei gh ts ( which we wi l l take the inverse of ) , from which γ pa ra mete rs ca n be es tim ate d using the Ande rse n and Gi l l ( 1982 ) model: In simulation studies, we as s es s ed both the log i stic r egr ession and the proportional rate model.Then, the FIPTM estimator s o lves the following equations: whe re E n s ta nds for the e mpirical mea n.How ev er, th at estim ator r equir es both the tr ea tme n t a nd the o bs erv ation models to be c orre ctly spe cifie d, which is not easy in practice.

N ovel estim ator
We propose the augme n ted AAIIW es tima tor ( which acr onym s ta nds for doub l y a ugmen ted an d do u bly in verse wei ght ed ) that is more flexible and allows two out of four different models to be misspec i fied while the estima tor r emain s con sis te n t.The es timator is dev elope d by finding the influenc e curv e of the es tima nd introduc e d in Se ct ion 2.4 ( Hines et al ., 2022 ) .The nove l e s tim ator is obtaine d by s o lving the fo llowing augme n te d v ersions of ( 4 ) : whe re the n uisa nce te rms ca n, for ins ta nce, be es tim ate d using pa ra metric models, with The la t te r model a rises whe n t aking the expect ation 5) .For the novel estimator, if using the proportional rate model for visits, then the baseline rate λ 0 (t ) in ( 2 ) m us t be es tim ate d before calculating the IIV wei gh ts, which was not the case with the FI PTM est imator.The IIV wei gh ts in the equations for the AAIIW are the inverse of E [dN i (t) | V i (t) ; γ] = λ 0 (t ) exp { γ T V i (t ) } .One can use the Breslow's estimator ( Cox, 1972 ) ( which we use in our sim ulation s tudies ) .
Table 2 shows the c ombin ations of c orre ctly spe c i fied models lea din g to a consis te n t AAIIW es timator.At leas t one of the two mode ls re lat ed t o confounde rs a nd one of the two models relat ed t o the o bs erv a tion pr ed ictor s m us t be c orre ctly spe c i fied.The es tima nd has a n efficie n t influe nce function if it is pathwise diffe re n t iable, i .e., if the univariate submodels are smooth in the pa ra mete rs, for the pos tul ated models ( s e e e g, Hines et al., 2022 ) .We herein ass ume th at the causal m argin al effe ct of treatme n t is pa th wise differ enti ab le.The deriv at ion of the est imator using the theory of influence functions and a proof of multiple ro bustnes s are in Web Appendices A and B , respe ctiv ely.The link betw e en the the ory of influenc e functions and the the ory on model -assis ted es t imat ion for our proposed estimator is in Web Appendix C .The c orre ct spe cification for pa ra metric model s i s elaborated in Web Appendix D .
We show in Web Appendix E that the AAI IW asymptot ic va ria nce, de rived as the v ari ance of its influence function ( Tsiatis, 2006 ) , is sm aller th an th at of the FIPTM when all n uisa nce models are c orre ctly spe c i fie d for both estim ators .In practic e, the AAIIW can be obtained by est imat ing the nuisance models para mete rs using e.g., log i stic r egr essions for the tr ea tme n t a nd visit models or a proportional rate model for the visits with coxph TA BLE 3 Est ima tors compar ed in simulation studies.

Estimator
OLS is a n ordina ry leas t squa re s e st imator, I PT is an I PT-wei gh ted es timator, DW is the ( doubly-wei gh ted ) FIPTM es tima tor fr om Coulombe et al . ( 2021 ) and AAI IW is the novel propose d estim ator.A means c orre ctly spe c i fied a nd a † symbol mea n s it is us ed as a nuis ance model in the estimator but it is wrongly spec i fied.
in R and two linear models for the outcome conditional on A i (t ) and K i (t ) or V i (t ) .The e stimate s can be plugged into the est imat ing equat ions of the AAIIW.Root s o lv ers ( s uch as uniroot in R ) can be used to estimate β 0 and that estimate be plugged into the se c ond e quation, s o lv e d for β 1 .
In Web Appendix F , we propose to relax the assumption that cen s oring occ ur s at random and to us e invers e pro bability of censoring wei gh ts ( IPCW ) ( Ro bin s a nd Finkels t ein, 2000 ) t o addr ess informa tive cen s oring.We further outline a multiply robust a pproach conside ring ce n s oring pred ictor s.The sce na rio with inform ativ e c en s oring i s al s o as s es s ed in simul ation s.

SIMUL ATION STUDY
In sim ulation s tudies, w e c ompare d four diffe re n t es timators detailed in Table 3 .
The DGM was str ongly inspir ed by simila r sim ulation s tudie s pre se n ted in B ůžková a nd Lumley ( 2009 ) , Coulombe et al. ( 2021Coulombe et al. ( , 2022 ) ) and is detailed more thoroughly in Web Appen dix G .The DGM included a s e t of confounders at baseline repea ted thr ough follo w -up, a time-va rying bina ry tr ea tme n t, a s e t of o bs erv a tion pr ed ictor s th at varie d in time, and irre gul ar o bs ervation of the outcome.The causal effect of tr ea tme n t was cons ta n t, i.e., w e c orre ctly spe cifie d the MSM in our simulations.The m ain res ults for 1000 simul ation s using a nonhomo geneous Pois s on rate to simulate the o bs erv at ion t imes of the outcome and a s amp le of size 1000 are prese n ted in the following Section 3.1 .In a nothe r s e t of simul ation s, we rep l ac e d the nonhomo geneous Pois s on rate with a nonhomo ge neous Be rnoulli probability for the observation indicator and used a log i stic regres sion in stead of the Ande rse n a nd Gi l l model to fit the probability of o bs erv a tion a t each time point.These r esults ar e pr ese n ted in Web Appendix H ( Web Figures 2 and 3 ) , along with the results under a s amp le of size 250 instead of 1000 ( Web Fig ure 1 ) and all Monte Carlo bi as es and mean squar e err ors ( We b Table 1 ) .In both s e ttin gs usin g either the Pois s on rate of the Bernoul li probabi lity, we t est ed four diffe re n t s e ts of γ pa ra meters in the o bs erv at ion model , including one s e t of zeros ( which we call "s e t 1" in the r esults ) , corr esponding to uninform ativ e obs erv ation .In another s en sitivity an alysis, w e as s es s ed the perform anc e of the propose d estim ator under inform ativ e c en s oring that depends on the visit pred ictor s V i (t ) .We compared IPCwei gh ted a nd more n aiv e estim ators th at do not addres s cen s oring.The simul ation s e tup is des cribed in Web Appe ndix G a nd the results from that a nalysis ( e mp irical b ias a nd mea n squa red e rror ) a re shown in Web Table 2 ( Web Appe ndix H ) a nd briefly d isc ussed in Section 3.1 .

Results
The distributions of 1000 e stimate s obt ained with each e stimator using a s amp le of size 1000 patients ar e pr esented in Figure 2 .A thorough d isc ussion of the results is given in Web Appendix H with more details on the pe rforma nce of each of the more naive estim ators .
The results a re ge ne rally as expe cte d .The AAI IW est imator is emp irically unb i as ed in all s c en arios ( 1 ) -( 4) for the o bs erv ation proc ess, whenev er using one of the four c ombin ations of corre ctly spe c i fied models show n in Table 2 or w he n all four models are c orre ctly spe c i fied.It exhibits small v ari ance when the tw o c ondition al outc ome mea n models a r e corr ectly spec i fied ( sc en ario b from Table 2 ) or, as expe cte d when all four models ar e corr ectly spec i fied.Results for the second s e t of simul ation s using the Bernoulli probability to simulate the o bs erv ation s, and those for a sample of size 250 are in Web Appendix H .As expe cte d , the est ima tors ar e mor e v ari ab le when using a s amp le size of 250, although the same pa t te rns in the compa rison of es timators are o bs erv e d ( Web Figure 1 ) .Similar results are o bs erv e d when using the Bernoulli probability instead of the Poisson rate for the simulation of observa tion indica tors ( Web Figur es 3 and 3 ) .The simul ation s using the Bernoul li probabi l ity d id not require the use of Breslow's estimator for the baseline rate, which may partly exp l ain the smaller v ari ances o bs erv e d ov erall ( e g, compa re Fi gure 1 and Web Figure 2 ) .
Results for the DGM with informative cen s oring that depended on the pred ictor s of visit were also as expected ( Web Table 2, Web Appendix H ) .The inform ativ e c en s oring did affect the emp irical b i as in s ome of the s c en arios t est ed.Adjustment via IPCW brought the es timates close r to the true causal effect, with a maximum bias that we n t from 0.31 to 0.14 after adjustme n t , for the AAII W estimator.The AAII W estim ator c ouple d with IPCW pe rformed pa rticula rly well whe n the tw o outc ome c ondition al mean models w ere c orre ctly spe c i fied, or when the  Acronyms: CI, c onfidenc e interval; IPT, inv ers e pro bability of tr ea tme n t; AIPW, augme n te d inv ers e pro bability of tr ea tme n t wei gh ted; IIV, inve rse in te nsit y of visit ; FIPTM, the flexib le invers e pro bability of tr ea tme n t a nd monitoring; AAIIW, the doubly augme n te d, doubly inv erse w ei gh ted.φ .Note we do not know the true data ge ne rating me ch anism for the tr ea tme n t me ch a nism in the a pp lication .† .This estimator us es a c orre ctly spe cifie d ge ne rating me ch anism for outc ome mis singnes s. ‡ .This estimator uses a wrongly spec i fied ge ne rating me ch anism for outc ome mis singnes s.
outcome model conditional on the confounders and the IIV wei gh ts we r e corr ectly spec i fied ( bias smalle r tha n 0.02 in all sc en arios ) .

M OTIVATIN G E X A MPLE
We applied the proposed AAI IW est imator a nd diffe re n t more n aiv e c omparat ors t o longitudinal da ta fr om the Ad d He alth study in the Unit ed Stat es ( Ha rris a nd Udry, 2022 ) .More details on that study and the analysis are av ail ab le in Web Appen dix I .We h av e ac c es s to d a ta fr om the first four waves of the Ad d He alth study, corre sponding to the years 1994-1995, 1996, 2001-2002, and 2008-2009, respe ctiv ely.Various types of information, including de mogra phics a nd health s tatus va riables we re c olle cte d via que stionnaire s fil led by the American ado les cents in this study.Our goal was to estimate the marginal causal effect of coun s eling on alco ho l con sumption bas ed on the question In the past y e ar, hav e y ou receiv e d psy c ho logica l or emotion a l co unseling? .The ass ume d D GM is shown in Fig ure 1 .Two challenges we wa n t ed t o consider in the analysis are the irregular observation of the outcome and, because the study is observat ional , the pote n tial confounding of the psy chother apy-alco ho l con sumption r ela tionship.We sele cte d sev e ral pote n tial confounde rs for that r ela tionship ( Web Appe ndix I ) .The a nalysis d atas e t c ontaine d s everal mis sing v alues.We us ed multip le imputation s by ch aine d equations ( Rubin, 1988 ) five times, to impute mis sing v alues in cov ari at es.The out come was defined using the question T h i n k of all the times you had a dri n k duri ng the past 12 months.How many dri n ks d i d yo u usua lly h av e e ach ti m e? .It con sists of a s elfas s es s ed n umbe r of drinks the ado les ce n t would consume, on averag e, e ach time they con sumed alco ho l, ran gin g from 0 to 90.In this application, the outcome was assessed at each of the four waves for everyone ( ie, it contained no missing value ) .To ass es s the adv a n tage of our approach, we simul ated mis singnes s in the outcome and as s es s ed the diffe re n t es timators in that s e tting, knowing the true underlying mis singnes s me ch anism .As suming that all pote n tial confounde rs as well as the mediator ( de pre ssive mood ) and the exposure ( coun s e ling ) a ffect the chance of obs erving the alco ho l con sumption out come, the out come o bs ervation ( ie, the opposite of missingness ) was simulated using a pre-spec i fied, inve n ted model ( W eb Appendix I ) .W e conducted the analysis usin g ea ch of the five imputed d atas e ts one by one.
We used Rubin's rule ( Rubin, 1976 ) to combine the final estima tes fr om all the estima tors compar ed, and 500 bootstrap s amp les to obtain c onfidenc e intervals ( CI ) .We fit a propensity score model and two diffe re n t proportional rate models for the o bs erv ation of the outc ome, one c orre ctly spe cifie d and one that was not c orre ctly spe c i fied ( as a function of the sinus of age and de pre ssive mood only ) .
An ordina ry leas t squa re s e s timator, a n IIV-wei gh ted es timator th at ac c ounts for the o bs erv a tion pr oc ess ( w e t est e d the tw o s e ts of the IIV wei gh ts ) , a doubly-wei gh ted es tima tor corr esponding to the FI PTM est im ator ( inc orporat ing the I PT wei gh ts based on our assumptions on the pote n tial confounde rs, a nd IIV wei gh ts-we tes te d the tw o s e ts of IIV wei gh ts ) , a nd the AAIIW estimator in which we incorporated the IPT wei gh ts a nd the two diffe re n t s e ts of the IIV wei gh ts we re compa red.We also added a comp le te d ata analysis in which an OLS, a n IPT-wei gh ted a nd a n augme n te d inv ers e pro bability of tr ea tme n t wei gh ted ( AIPW ) estim ators w ere c ompute d on the d atas e t with no mis sing d ata for the outcome.Some diffe re nc es w ere found across the tw o expos ure groups in the first imputed d atas e t, which indicates pote n tial confounding ( Web Appendix J, Web Table 3 ) .In the outcome observation model, we also found modest differences in female sex and smoking status betw e en those for whom the alcohol c ons umption was o bs erv e d and the others ( Web Appendix J , Web Table 4 ) .After IIV wei gh ting, mos t diffe re nces va nished ( Web Table 4 ) .
Both the adjus tme n t for confounding a nd the one for outcome mis singnes s bring the e stimate s for the m argin al effe ct of expos ure to c oun s eling tow a rd the n ull.The es tim ator th at le d to the close st e stimate s to the comp le te d ata a nalysis ( poin t es timate 0.35 with the AIPW, Table 4 ) is the AAI IW est im ator, which le d to poin t es timates of 0.40 and 0.39 when using the c orre ct or the wrong IIV wei gh ts, respe ctiv ely.The FI PTM est im ator le d to point e stimate s of 0.36 and 0.72, respe ctiv ely ( Table 4 ) , with the estimator using the wrong IIV wei gh ts leading to the estimate further away from the gold s ta nda rd poin t es tima te.Our r esults indica te tha t in a s e tting in which we would not know the true o bs erv ation me ch anis m, the AAIIW e stimator might sti l l lead to a n es timate of the causal effect closer to the comp le te d ata a nalysis, while the FIPTM is more at risk of being bi as ed if its inv erse w eights ar e wr ongly spec i fied.Our pr oposed appr oach allows a djustin g for previous ( o bs erv e d ) treatme n ts or outcomes as pote n tial confounde r s or visit pred ictor s, but it cannot address s e ttings in which a previous outcome ( that is not o bs erved ) affects the o bs erv a tion of any futur e outcome.In this application, we did not include previous outcomes in the adjus tme n t s e t for tha t r eason, even if pr evious outc ome values w ere av ail ab le.The proposed AAI IW est imator can be used to estimate a time-fixed averag e tre atme n t effect.In this application, it is pos sib le that the causal effect of the ra p y on alco ho l con sumption changes in time, with e.g., a gr ea te r be nefit of the ra p y at the beginning of follo wup, but we e s timated a n "ave raged" ove r all times tr ea tme n t effect, β 1 .The true tr ea tme n t effect could vary in time.A lengthier d isc ussion on the study results is given in Web Appendix J .

DIS CUSS ION
This work proposed the firs t m ultip ly ro bus t es timator for the causal m argin al effe ct of tr ea tme n t a ddressin g confoundin g and irre gular visits, th at is consis te n t whe n only two out of four n uisance models, one r ela t ed t o confounde rs a nd one to visit pred ictor s, ar e corr ectly spec i fied .In addit ion to being more robust than the FIPTM, the AAIIW estimator i s al so the most efficient estimator in its s emiparame tric cl as s.In simul ation studies, it w as demon strat ed t o be robus t a nd e mpirically as efficie n t as the FIPTM when the two wei gh t models a r e corr ectly spec i fied but it could be much more efficient in some other sc en arios .
In an application to the Ad d He alth study in the United St ate s, we found a diffe re nce betwee n more naive es timators a nd the multip ly ro bust AAI IW est imator in the est imat ion of the causal m argin al effe ct of the ra p y coun s eling on alco ho l con sumption, and the proposed estim ator le d to the e stimate s th at w ere the closest to a gold s ta nda rd found with the complete dataset.It is pos sib le, how ev er, th at unmeas ure d c onfounding rem ains .Sensitivity analyses can be used to assess the effect of unmeasured confounding or visit pred ictor s tha t wer e not accounted properly in the estimator ( see eg, McCulloch and Neuhaus 2020 for diagnostics on visit irregularity when visit times may depend on the outcome values, or Va nde rWeele a nd Arah 2011 for se nsitivity analyses that address unmeas ure d c onfounding ) .
The consis te ncy of our proposed estima tor r e lie s on spec i fic c ombin ations of c orre ctly spe c i fied n uisa nce models listed in Table 2 and some classical causal assumptions me n tioned in Section 2 , includ ing cond ition al exch ang e a bility.S ee Web Appendi x K for some re c ommendations on the ide n t ificat ion of adjustme n t s e ts.The propos ed approach als o re lie s on the ass ume d MS M. We as s ume in this w ork th at the outc ome is r ela t ed t o the tr ea tme n t at time t by a constant parameter ( causal effect ) β 1 .Thus, our working model, the ass ume d MSM, is only c orre ctly spec i fied i f the tr ea tme n t causal effe ct is c ons ta n t, i.e., if it is the sa me for a ny time t.If it is not, the n the es tim ate d effe ct c orresponds to the closest time-fixed effect to the true , time-varying causal effect and acts as a summary of the true causal r ela tionship if all n uisa nce models ar e corr ectly spec i fied ( Ne u geb a ue r a nd va n de r Laa n, 2007 ) .Furthe r more, i f the working MSM model is not c orre ctly spe cifie d, causal in te rpr eta tion is more d iffic ult, as the estim ate d effe ct is av erage d ov e r all time poin ts a nd does not r epr ese n t the causal effect of tr ea tme n t at time t.That effect ca n ins tead be in te rpreted as the averag e tre atme n t effect ove r the e n tire follo w -up period, if one follow e d a c ons ta n t tr ea tme n t course ( A i (t ) = 1 for all t, or A i (t ) = 0 for all t) , but it be c omes h a rde r to in te rpre t if one fo llo ws a tre atme n t course with tr ea tme n t switches .In s uch s e ttings, nonpa ra me tric MS M s uch as propose d in Ne u geb a ue r a nd va n de r Laa n ( 2007 ) could be preferable to estimate causal curves as a function of time, or the tr ea tme n t a nd the visit proces s es could be modeled jointly to ackno wledg e the lack of g e ne rali zab ility of the effect at one time, to other times when there is no visit ( see eg, Ro bin s e t al. 2008 ; Ne u geb a uer et al. 2017 , who d isc ussed ide n t ificat ion of optimal tr ea tme n t a nd visit s trategies unde r join t models for the tw o proc es s es ) .

FIG URE 1
FIG URE 1 Caus al di agram i l lustrating the ass ume d DGM. ( A ) Ass ume d caus al di agram a t time t in pa tie n t i , pos tulat ed t o be common across all patie n ts. ( B ) Caus al di agra m for the Add Health s tudy da ta a t time t common acr os s all ado les ce n ts.
2. ( a ) positivity of tr ea tme n t, mea ning that a nyone should h av e a ch anc e of re c eiving any of the tw o tr ea tme n t opt ions, and ( b ) posit ivity of observation, s uch th at patie n ts had a chance to have their outcome o bs erved at any time given their cha racte ris tics.3. Condition al exch ang e ability, with includes ( a ) no unme as ure d c onfounder in the o bs erv e d s e t K i (t ) ; and ( b ) indepe nde nce of the o bs erv a tion indica tors with other variables in the analysis conditional on the visit pred ictor s V i (t ) .

FIGUR E 2
FIGUR E 2 R esults of the sim ulation s tudies with a sa mple size of 1000 using a nonhomogeneous Poisson rate to simulate the o bs erv ation ind icator s and the Andersen and Gi l l model with Breslow estimator to estimate the IIV wei gh ts.Each boxplot r epr ese n ts the distribution of 1000 e stimate s for the corre sponding e stim ator.The dashe d line r epr ese n ts the gold s ta nda rd , i .e., the true value for the m argin al effe ct of expos ure th at e quals to 1. Diffe re n t s tre ngth s of the visit proces s on cov ari a tes ar e r epr ese n ted with sc en arios ( A ) γ = (0 , 0 , 0 , 0 , 0 , −5) ( ie, no bias due to the visit process expe cte d ) ; ( B ) γ = (0 . 5 , 0 .3 , −0 .5 , −2 , 0 , −3) ; ( C ) γ = (0 . 5 , −0 .5 , −0 . 2 , −1 , 1 , −3) ; and ( D ) γ = (−1 , −0 .8 , 0 . 1 , 0 .3 , −1 , −3) .OLS, ordinary least squares; IPT, inverse probability of tr ea tme n t wei gh ts; a nd DW, doubly-wei gh ted estimator which corresponds to the FIPTM from Coulombe et al. ( 2021 ) ; AAIIW: The novel doubly augme n ted, doubly wei gh ted es timator.The subscripts c , nc , iptc , and iivc , respe ctiv ely mean all c orre ct, all not c orre ct, only IPT c orre ct, and only IIV c orre ct in the n uisa nc e models .The subscripts s.a to s.d refer to sc en arios ( A ) -( D ) in Table 2 of the ma n uscript.
is used to define our es tima nd.De note b y Y 1 i (t ) the pote n tial outcomes of individual i at time t if they

TABLE 1
Causal assumptions r equir ed for the proposed estimator to be consis te n t.
d t the martingale residual for the o bs erv a tion pr oc ess .The c ondition al outc ome mean models in the augme n t ed t erms are μ a

TABLE 4
Comp le t e out come data ( t op ) a nd irregula rly obse rv e d outc ome da ta ( bot tom ) estima tes ( 95% boots tra p pe rce n tiles CI ) of the m argin al effe ct of c oun s eling on the ave rage n umbe r of alco ho lic bev erages c ons ume d, Ad d He alth study, UnitedSt ate s, 1996St ate s,  -2008.   .